//ionization
#ifndef NONLINEARRESPONSE_H 
#define NONLINEARRESPONSE_H 
#include <cmath>
#include <complex>
#include <iostream>
#include <cstdio>
#define complex std::complex<double>

#include "./Constants.h"

using std::cout;
using std::endl;

class NonlinearResponse {
private:
    int m_nt, m_dt;
    double eta2 ; //optical Kerr coefficient in air (a.u.)    
    double Ng ;       // density_of_neturals: 2.4*10^19/cm^3  = 3.55609e-6 a.u. 
    double C = Constants::c;
    double eps0 = Constants::epsilon_0;

    //Ip: First-order ionization potential of nitrogen molecules, 
    //U_H: ionization potential of hydrogen atoms, 
    double Ip = 0.5;     //N_2氮气分子一阶电离势 15.581/27.2114=0.5726
    double U_H = 0.5;    //H氢原子电离势s
    double U = Ip / U_H;

    double Eabs(double E);
    double ADK (double E_abs);
    double current_Jabs (double Wt, double Net, double E_abs, double E_field);
    double current_Jp(double Net0, double E0, double Net1, double E1);
	double polarization(double E_field);
   
public:
    NonlinearResponse(int nt, double dt, double kerr, double initialDensity) :
             m_nt(nt), m_dt(dt), eta2(kerr), Ng(initialDensity){
    }

    virtual ~NonlinearResponse() {  
    }
 
    void operator()(double* ionization_rate, 
                    double* electron_density,
                    double* nonlinear_J,
                    double* nonlinear_P,
                    double* electric_field);
    
};

double NonlinearResponse::Eabs(double E){
    return fabs(E) + 1e-17;       //在fabs(Ereal) + a,以防E=0时不收敛
}
double NonlinearResponse::ADK (double E_abs){
    return 4 * (pow(U, 2.5)/E_abs) * exp(-(2 * pow(U, 1.5))/(3 * E_abs));
}

double NonlinearResponse::current_Jabs (double Wt, double Net, double E_abs, double E_field){
    double Jabs_value = (Wt * Net * Ip * E_field) / (std::pow(E_abs, 2));
    return Jabs_value;
}

double NonlinearResponse::current_Jp(double Net0, double E0, double Net1, double E1){
    double Jp_value = (Net0 * E0 + Net1 * E1) * m_dt * 0.5;
    return Jp_value;
}

double NonlinearResponse::polarization(double E_field){
    double P3_value = (eps0 * eta2 * std::pow(E_field, 3) / C);
    return P3_value;
}

void NonlinearResponse::operator()(double* ionization_rate, 
                                   double* electron_density,
                                   double* nonlinear_J,
                                   double* nonlinear_P,
                                   double* electric_field){
    
    double integral_gama = 0.;
    double Jp =0;
    double Jabs =0;
    for(int num_t = 0; num_t < m_nt; num_t++){
        double E = electric_field[num_t];
        double E_1 = electric_field[num_t-1];
        double E_abs = Eabs(E); 
        double Wt = ADK(E_abs);
        ionization_rate[num_t] = Wt;
        double Wt_1 = ionization_rate[num_t-1];

        integral_gama += (num_t == 0)? 0.: ((Wt + Wt_1) * m_dt * 0.5); 
        
        double Net   = (num_t == 0)? (Ng * Wt ): (Ng *  (1. - std::exp(- integral_gama)));
        electron_density[num_t] = Net;
        double Net_1 = electron_density[num_t-1];
        
        Jp += (num_t == 0) ? 0 : current_Jp(Net_1, E_1, Net , E);
        Jabs = current_Jabs(Wt, Net, E_abs, E);
        
        nonlinear_J[num_t] = (Jabs + Jp) * (4. * M_PI  / (std::pow(C, 2)));
        nonlinear_P[num_t] = polarization(E);
    }
}
#endif
